rm(list=ls())
#setwd("/Main Text-Figures/Figure4/")

########################################

alpha <- 1
k <- 0.2
N <- 100
v <- 500
nleg <- 50
nsim <- 100

diff.left.cohesion.sim <- NULL
diff.right.cohesion.sim <- NULL
cohesion.right.all.sim <- NULL
cohesion.right.rcv.sim <- NULL
cohesion.left.all.sim <- NULL
cohesion.left.rcv.sim <- NULL
sim.info.all <- NULL

for(z in 1:nleg){

phi <- runif(1,.5,.8)
  npL <- round((1-phi)*N,digits=0)
  npR <- round((phi)*N,digits=0)

#heterogeneity scores
d1 <- runif(n=1,min=0.5,max=2)
d2 <- d1

########################################

xL1 <- c(0,runif(npL-1,min=0-d1,max=d1)) #Going to use this that i've put the party leaders in the first and last spots
xL2 <- xL1#c(0,runif(npL-1,min=0-d2[k],max=d2[k]))
xR1 <- c(runif(npR-1,min=1-d1,max=1+d1),1) #Going to use this that i've put the party leaders in the first and last spots
xR2 <- xR1#c(runif(npR-1,min=1-d2[k],max=1+d2[k]),1)

left.party <- as.data.frame(cbind(xL1,xL2,rep(0,npL)))
colnames(left.party) <- c("x1","x2","partyid")
right.party <- as.data.frame(cbind(xR1,xR2,rep(1,npR)))
colnames(right.party) <- c("x1","x2","partyid")
all.leg <- as.data.frame(rbind(left.party,right.party))

for(j in 1:nsim){
  b1 <- runif(v,0,1)
  b2 <- b1
  sq1 <- runif(v,0,1)
  sq2 <- sq1#runif(v,0,1)
  
  hypo.pro.1d <- cbind(rep(1,length(b1)),rep(d1,length(b1)),b1,sq1)
  colnames(hypo.pro.1d) <- c("dim","d","b","sq")
  hypo.pro.2d <- cbind(rep(2,length(b2)),rep(d2,length(b2)),b2,sq2)
  colnames(hypo.pro.2d) <- c("dim","d","b","sq")
  hypo.pro.all <- as.data.frame(rbind(hypo.pro.1d,hypo.pro.2d))
  hypo.pro.all$id <- seq(1,nrow(hypo.pro.all),1)

########################################
#votes?
    vote.func1 <- function(x1,b1,sq1){
      as.numeric(abs(x1-b1)  < abs(x1-sq1) )
    }
    vote.func2 <- function(x2,b2,sq2){
      as.numeric(abs(x2-b2)  < abs(x2-sq2) )
    }

    out <- matrix(data=NA,nrow=N,ncol=v*2);dim(out)
    for(i in 1:N){
      for(j in 1:v){
        out[i,j] <- vote.func1(all.leg$x1[i],b1[j],sq1[j])
      }
    }
    for(i in 1:N){
      for(j in (v+1):(2*v)){
        out[i,j] <- vote.func2(all.leg$x2[i],b2[j-v],sq2[j-v])
      }
    }
########################################
#pass?
    
    passes <- as.numeric(colSums(out) > .5*N)
    hypo.pro.all <- as.data.frame(cbind(hypo.pro.all,passes))
    for (i in 1:nrow(hypo.pro.all)){
      hypo.pro.all$outcome[i] <- if(hypo.pro.all$passes[i]==0) hypo.pro.all$sq[i] else hypo.pro.all$b[i]
    }

########################################
#unity score?
    left.aligned <- rep(NA,length=ncol(out))
    for(i in 1:ncol(out)){
      left.aligned[i] <- round((sum(as.numeric(out[2:npL,i]==out[1,i])))/npL,4)
    }
    right.aligned <- rep(NA,length=ncol(out))
    for(i in 1:ncol(out)){
      right.aligned[i] <- round((sum(as.numeric(out[(npL+1):(N-1),i]==out[N,i])))/npR,4)
    }

    out.cohesion <- matrix(data=NA,nrow=N,ncol=v*2)
    for(j in 1:(v*2)){
      for(i in 1:npL){
        out.cohesion[i,j] <- if (out[i,j] == out[1,j]) {
          left.aligned[j] - right.aligned[j]
        } else {
          0
        }
      }
    }
    for(j in 1:(v*2)){
      for(i in (npL+1):N){
        out.cohesion[i,j] <- if (out[i,j] == out[N,j]) {
          right.aligned[j]-left.aligned[j]
        } else {
          0
        }
      }
    }
    
    out.cohesion.abs <- matrix(data=NA,nrow=N,ncol=v*2)
    for(j in 1:(v*2)){
      for(i in 1:npL){
        out.cohesion.abs[i,j] <- if (out[i,j] == out[1,j]) {
          left.aligned[j] 
        } else {
          0
        }
      }
    }
    for(j in 1:(v*2)){
      for(i in (npL+1):N){
        out.cohesion.abs[i,j] <- if (out[i,j] == out[N,j]) {
          right.aligned[j]
        } else {
          0
        }
      }
    }

########################################
#RCV requested?
    
    out.cohesion.kmat <- matrix(data=NA,nrow=N,ncol=v*2)
    for(i in 1:N){
      for(j in 1:(v*2)){
        out.cohesion.kmat[i,j] <- as.numeric(out.cohesion[i,j] > k)
      }
    }
    hypo.pro.all$left.cscore <- out.cohesion[1,] 
    hypo.pro.all$right.cscore <- out.cohesion[N,]
    
    hypo.pro.all$left.cscore.abs <- out.cohesion.abs[1,] 
    hypo.pro.all$right.cscore.abs <- out.cohesion.abs[N,]
    
    hypo.pro.all$left.rcv.request <- out.cohesion.kmat[1,]
    hypo.pro.all$right.rcv.request <- out.cohesion.kmat[N,]
    
    hypo.pro.all$rcv.request <- as.numeric(colSums(out.cohesion.kmat)>0)

########################################
#proposed?

    for(i in 1:nrow(hypo.pro.all)){
      hypo.pro.all$proposing.party[i] <- if(hypo.pro.all$b[i] < hypo.pro.all$sq[i]) 0 else 1 
    }
    hypo.pro.all$proposed.left <- rep(0,length(nrow(hypo.pro.all)))
    for(i in 1:nrow(hypo.pro.all)){
      if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i]-k > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i]) > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$b[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==0 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.left[i] <- as.numeric(alpha*-abs(0-hypo.pro.all$sq[i])+hypo.pro.all$left.cscore[i] > alpha*-abs(0-hypo.pro.all$sq[i]))
      } else {
        0
      }
    }
    hypo.pro.all$proposed.right <- rep(0,length(nrow(hypo.pro.all)))
    for(i in 1:nrow(hypo.pro.all)){
      if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$right.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i]-k > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$rcv.request[i]==0 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i]) > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==1) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$b[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else if (hypo.pro.all$proposing.party[i]==1 & hypo.pro.all$left.rcv.request[i]==1 & hypo.pro.all$passes[i]==0) {
        hypo.pro.all$proposed.right[i] <- as.numeric(alpha*-abs(1-hypo.pro.all$sq[i])+hypo.pro.all$right.cscore[i] > alpha*-abs(1-hypo.pro.all$sq[i]))
      } else {
        0
      }
    }

    hypo.pro.all$proposed <- hypo.pro.all$proposed.left+hypo.pro.all$proposed.right

    all.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1),]
    all.proposals <- all.proposals[which(all.proposals$rcv.request==0),]
      all.proposals.left <- all.proposals[which(all.proposals$proposed.left==1),]
      all.proposals.right <- all.proposals[which(all.proposals$proposed.right==1),]
    rcv.proposals <- hypo.pro.all[which(hypo.pro.all$proposed==1 & hypo.pro.all$rcv.request==1),]
      rcv.proposals.left <- rcv.proposals[which(rcv.proposals$proposed.left==1),]
      rcv.proposals.right <- rcv.proposals[which(rcv.proposals$proposed.right==1),]

########################################

    vote.matrix.all <- out[,all.proposals$id]
      vote.matrix.all.left <- out[,all.proposals.left$id]
      vote.matrix.all.right <- out[,all.proposals.right$id]
    vote.matrix.rcv <- out[,rcv.proposals$id]
      vote.matrix.rcv.left <- out[,rcv.proposals.left$id]
      vote.matrix.rcv.right <- out[,rcv.proposals.right$id]
    
########################################
#Cohesion
    
    cpd.all <- as.data.frame(c(all.proposals$left.cscore.abs,all.proposals$right.cscore.abs))
    cpd.all$id <- c(rep("L",length(all.proposals$left.cscore.abs)),rep("R",length(all.proposals$right.cscore.abs)))
    colnames(cpd.all) <- c("cscore","id")
    
    cpd.rcv <- as.data.frame(c(rcv.proposals$left.cscore.abs,rcv.proposals$right.cscore.abs))
    cpd.rcv$id <- c(rep("L",length(rcv.proposals$left.cscore.abs)),rep("R",length(rcv.proposals$right.cscore.abs)))
    colnames(cpd.rcv) <- c("cscore","id")
    
    cpd.right <- as.data.frame(c(all.proposals$right.cscore.abs,rcv.proposals$right.cscore.abs))
    cpd.right$id <- c(rep("ALL",length(all.proposals$right.cscore.abs)),rep("RCV",length(rcv.proposals$right.cscore.abs)))
    colnames(cpd.right) <- c("cscore","id")

    cpd.left <- as.data.frame(c(all.proposals$left.cscore.abs,rcv.proposals$left.cscore.abs))
    cpd.left$id <- c(rep("ALL",length(all.proposals$left.cscore.abs)),rep("RCV",length(rcv.proposals$left.cscore.abs)))
    colnames(cpd.left) <- c("cscore","id")

    diff.left.cohesion <- mean(all.proposals$left.cscore.abs,na.rm=TRUE) - mean(rcv.proposals$left.cscore.abs,na.rm=TRUE)
    diff.right.cohesion <- mean(all.proposals$right.cscore.abs,na.rm=TRUE) - mean(rcv.proposals$right.cscore.abs,na.rm=TRUE)
    cohesion.right.all <- mean(all.proposals$right.cscore.abs,na.rm=TRUE)
    cohesion.right.rcv <- mean(rcv.proposals$right.cscore.abs,na.rm=TRUE)
    cohesion.left.all <- mean(all.proposals$left.cscore.abs,na.rm=TRUE)
    cohesion.left.rcv <- mean(rcv.proposals$left.cscore.abs,na.rm=TRUE)
  }

  diff.left.cohesion.sim <-  c(diff.left.cohesion.sim,diff.left.cohesion)
  diff.right.cohesion.sim <- c(diff.right.cohesion.sim,diff.right.cohesion)

  cohesion.right.all.sim <- c(cohesion.right.all.sim,cohesion.right.all)
  cohesion.right.rcv.sim <- c(cohesion.right.rcv.sim,cohesion.right.rcv)

  cohesion.left.all.sim <- c(cohesion.left.all.sim,cohesion.left.all)
  cohesion.left.rcv.sim <- c(cohesion.left.rcv.sim,cohesion.left.rcv)

  sim.info <- c(nrow(all.proposals),nrow(rcv.proposals),phi,d1)
  sim.info.all <- rbind(sim.info.all,sim.info)
}

sim.info.sum <- sim.info.all

abs.left <- abs(diff.left.cohesion.sim)
abs.right <- abs(diff.right.cohesion.sim)
abs.all <- c(abs.left,abs.right)
table(abs.all>.1)

all <- c(diff.left.cohesion.sim,diff.right.cohesion.sim)
table(all>0)

sim.info.left <- cbind(sim.info.all,diff.left.cohesion.sim)
colnames(sim.info.left) <- c("non.rcv.count","rcv.count","phi","d","diff")
sim.info.right <- cbind(sim.info.all,diff.right.cohesion.sim)
colnames(sim.info.right) <- c("non.rcv.count","rcv.count","phi","d","diff")
sim.info.all <- rbind(sim.info.left,sim.info.right)

sim.info.all <- as.data.frame(sim.info.all)
sim.info.all$abs.diff <- abs(sim.info.all$diff)

sim.info1 <- sim.info.all[which(sim.info.all$phi >=.35 & sim.info.all$phi <= .65 & sim.info.all$d >=.75 & sim.info.all$d <= 1.5),]
summary(sim.info1)
sim.info2 <- sim.info.all[which( (sim.info.all$phi< .35 | sim.info.all$phi > .65) & (sim.info.all$d < .75 | sim.info.all$d > 1.5)),]
summary(sim.info2)
sd(sim.info2$abs.diff)

load("cohesion-fig4-data.RData")

#pdf("Fig4_.pdf",width=6,height=6)
par(mfrow=c(1,1),mar=c(5,5,1,2))
hist(all,border="white",col="grey60",
     xlab="Non-RCV Cohesion - RCV Cohesion",
     main="",xaxt='n',
     xlim=c(-.25,.25),ylim=c(0,12),yaxt='n',
     breaks=seq(-.5,.5,.025),ylab="")
mtext(side=3,line=1,"",cex=1.5)
axis(side=2,at=seq(0,12,2),labels=seq(0,12,2),las=1)
axis(side=1,at=seq(-.25,.25,.05),labels=c("",-.2,"",-.1,"",0,"",.1,"",.2,""),las=1)
mtext(side=2,line=3,"Frequency")
#dev.off()


#pdf(".pdf",width=8,height=1)
#par(mar=c(0,0,0,0))
#plot(0,0,xlim=c(-1,1),ylim=c(.5,1),col="white",ylab="",xlab="",xaxt='n',yaxt='n',bty='n')
#legend(-1,1,legend=c("Full Sample","",expression(paste("Sample: d"%in%"[.75,1.5]")),
#                     expression(paste("               ",phi%in%"[.35,.65]"))),
#       pch=c(15,0,15,0),cex=2, col=c("black","white","grey60","white"),bty='n',ncol=2)     
#dev.off()


colnames(sim.info.sum) <- c("no.rcv","rcv","phi","d")
sim.info.sum <- as.data.frame(sim.info.sum)
summary(sim.info.sum)
sd(sim.info.sum$d)
